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We investigate the stability properties of breather sohton trains in a three-dimensional 
Bose-Einstein Condensate with Feshbach Resonance Management of the scattering 
length. This is done so as to generate both attractive and repulsive interaction. The 
condensate is confined only by a one dimensional optical lattice and we consider 
both strong, moderate, and weak confinement. By strong confinement we mean a 
situation in which a quasi two dimensional soliton is created. Moderate confinement 
admits a fully three dimensional soliton. Weak confinement allows individual soli- 
tons to interact. Stability properties are investigated by several theoretical methods 
such as a variational analysis, treatment of motion in effective potential wells, and 
collapse dynamics. Armed with all the information forthcoming from these meth- 
ods, we then undertake a numerical calculation. Our theoretical predictions are 
fully confirmed, perhaps to a higher degree than expected. We compare regions of 
stability in parameter space obtained from a fully 3D analysis with those from a 
quasi two-dimensional treatment, when the dynamics in one direction are frozen. 
We find that in the 3D case the stability region splits into two parts. However, 
as we tighten the confinement, one of the islands of stability moves toward higher 
frequencies and the lower frequency region becomes more and more like that for 
quasi 2D. We demonstrate these solutions in direct numerical simulations and, im- 
portantly, suggest a way of creating robust 3D solitons in experiments in a Bose 
Einstein Condensate in a one-dimensional lattice. 
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Feshbach resonance 

1. Introduction 

The creation of Bose-Einstein condensates (BEC) in vapours of alkali metals offers 
a splendid opportunity to investigate the wide range of nonlinear effects invloving 
atomic matter waves. A particular challenge is to develop methods for creating and 
controlling matter- wave solitons. Both dark and bright one dimensional solitons in 
harmonic traps have been observed (Burger et. al. 1999; Denschlag et. al. 2000; 
Strecker et. al. 2002; Khaykovich et. al. 2002; Eiermann et. al. 2004). A promising 
approach to obtaining multidimensional solitons consists in varying the scattering 
length (SL) of interatomic collisions. This can be achieved by means of sweeping an 
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external magnetic field through the zero-SL point close to the Feshbach resonance 
(Donley et. al. 2001; Inouye et. al. 1998; Fedichev et. al. 1996; Saito & Ueda 2002; 
Thcis et. al. 2004) . The application of an ac magnetic field may induce a periodic 
modulation of the SL, opening the way to so called 'Feshbach-resonance manage- 
ment' (FRM) (Kevrekidis et. al. 2003). A noteworthy FRM- induced effect is the 
possibility of creating self-trapped oscillating BEC solitons (breathers) without an 
external trap in the 2D case. The underlying idea is that fast modulations create 
an effective potential on a slower timescale. This potential can stabilize the soliton. 

The mathematical model for a BEC based on the Gross-Pitaevskii equation 
(GPE) with harmonic modulation of the SL was investigated by Saito & Ueda 
(2002, 2004), Abdullaev et. al. (2003), and Montesinos et. al. (2004). The conclu- 
sion was that FRM enables us to stabilize 2D breather solitons even without the use 
of an external trap. According to these references, 3D breathers conversely require 
at least a tight, one dimensional harmonic trap, practically reducing the problem 
to 2D (Gorlitz et. al. 2001). Later on wc will call this approach a quasi two dimen- 
sional treatment (Q2D). It will be defined more precisely below. Stabilization of 3D 
droplets is also possible by including dissipative effects (Saito & Ueda 2004). 

Trippenbach et. al. (2004) demonstrated that 3D solitons can be stabilized by 
a combination of FRM and a ID optical lattice (ID OL), instead of a ID harmonic 
trap. This issue has practical relevance, as a ID OL can easily be created by illu- 
minating the BEC by a pair of counterpropagating laser beams so that they form a 
periodic interference pattern (Stecher et. al. 1997). Incidently, it is easier to realize 
a tight confinement configuration in an optical lattice than in a harmonic trap. 
Hence this environment may well be more friendly for creating quasi 2D solitons. 
The lattice will be weak or strong depending on the intensity of the laser light. The 
combined OL-FRM stabilization of 3D solitons is possible even in a weak lattice, 
when atoms confined in different cells do interact. 

By analyzing the stability charts in configuration space we find two distinct 
regions where stable solutions exist. The first of these regions has its counterpart 
in the Q2D treatment. The other region appears when the frequency of modulation 
exceeds the lowest excitation frequency of the confining potential. It is not present 
in the Q2D treatment; and so can only correspond to fully 3D solitons. In the limit 
of tight confinement, the latter region moves up to extremely high frequencies and 
the Q2D stability chart is recovered. 

2. Theoretical approach 

We describe our system by the GPE in reduced units, including a time-dependent 
(FRM-controUed) nonlinear coefficient g{t) and an external potential U{r,t) 

i^t= [-(l/2)V^ + C/(r,t)+5(i)|^|2]^. (2.1) 

Initially the BEC is in the ground state of a radial (2D) parabolic trap with fre- 
quency Lu_i_, supplemented, in the longitudinal direction, by 'end caps', induced by 
transverse light sheets. The configuration is very similar to that used to create soli- 
ton trains in a Li'' condensate (Khaykovich et. al. 2002). A ID lattice potential in 
the axial direction is adiabatically turned on from £ = 0to£: = £/, see figure 1. 



Article submitted to Royal Society 



Stability analysis of three dimensional breather solitons... 



3 




/■ 



Figure 1. The time dependence of the nonhnear coefficient, g, switching-off function, f{t), 
and optical-lattice strength, e, in the numerical experiment. This experiment leads to 
the establishment of stable 3D breathing solitons supported by the combination of the 
quasi- ID lattice and Feshbach-resonance management (FRM). The shaded area indicates 
rapid oscillations of g{t), which account for the FRM. 

Thus, the full potential is, with period normalized to tt 



where g is the radial variable in the plane transverse to z, and the axial 'end-cap' 
potential, Uo{z) is approximated by a sufBciently deep one dimensional rectangular 
potential well. The width of the well determines the number of peaks in the final 
structure. The f{t) is a switch function (see figure 1). 
The nonlinear interaction coupling is described by 



Initially gi{0) = and g{0) = go{0) > 0. At some moment ti, we begin decreasing 
go{t) linearly. It vanishes at time t2, and remains zero up to ts, when we start to 
gradually switch on the rapid FRM modulation of the SL. In the interval [^3,^4], 
go{t) decreases linearly from zero to a negative value gof and the amplitude of the 
modulation gi{t) increases from zero to gif, see figure 1. Simultaneously, both the 
radial confinement and end-caps arc gradually switched off (sec the behaviour of the 
f{t) function in figure 1). At times t > ti, g{t) oscillates with a constant amplitude 
g\f around a negative average value 50/- Consequently, a soliton so created, if any, 
is supported by the combination of the ID lattice and FRM. Note that we choose 
a set of specific ramp functions go and gi and they grow rather rapidly in time. In 
general the stability diagram may depend on the shape and rapidity of the ramp 
functions. We did not investigate this issue in detail, but we believe that our results 
are to some extend universal. The reason we start from a positive value of go is 
so that then the atoms can be uniformly distributed in the cells. This will allow 
stabilization of solitons in all filled cells simultaneously. 

Numerical experiments following the path outlined in figure 1 indicate that it 
is possible to create stable 3D solitons (see inset to figure 2). Before displaying 
the results, we first resort to the variational approximation (VA) in order to pre- 
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diet conditions on the modulation frequency and the size of the negative average 
nonUnear coefficient gof, necessary to support 3D sohtons. 



3. Variational analysis 

The VA can be appHed to the description of BEC dynamics under diverse circum- 
stances (e. g. Saito & Ueda 2003; Montesinos et. al. 2004; Baizakov et. al. 2003). 
Equation (2.1) is derived from the Lagrangian density 

£ = i{^p:^ - rt^) - iv.i" - iV'.i' - 9{tm'' - 2t/|v^i^ (3.1) 

We use the VA for t > tj, and choose a hybrid ansatz composed of an hyperbolic 
secans function and a Gaussian for the solution for one lattice cell (calculations with 
two Gaussians give slightly inferior results, in terms of agreement with numerics). 
The amplitude is A{t)^ the overall phase is (j), the radial and axial widths are W{t) 
and V{t) respectively, and 6(t) and 0{t) are the corresponding chirps 

V(r,i) = Asech(e/W^)e[-'^^'-(i/2^'+''^)^'+H. (3.2) 

The reduced Lagrangian can be found upon substituting (3.2) into (3.1) and inte- 
grating over space. It is 



A^gh , eh f{t)LolW^h h ^ 



2V2 e^' 2 4t/2 



(3.3) 



where h = 27rln2, h = (9/4)7rC(3), and h = (7r/3)(41n2 - 1). By varying this 
reduced Lagrangian with respect to wc obtain the constant E = A^W^V = 
(27r^/2 ln2)~^ /^gjJVpdr = (2n7r^/2 ln2)~^, where the integral extends over one 
cell of the lattice, and n is the number of occupied lattice cells. Notice that the 
total number of atoms is included in the definition of the nonlinear coupling g{t), 
and the total wavefunction is normalized to unity. When the other four variational 
equations are derived, we can deduce two dynamical equations for the widths by 
substituting for the chirps 

b = W/W, (3.4) 
/? = V/V, (3.5) 



to obtain 



i3-4.,yexp(-F^) + |0, (3.7) 



V = 

where Ji = {h - h)/^, J2 = Eh/{y/2h), and J3 = Ehl{^h). k\\ J's are 
positive. These equations describe the dynamics of a single peak, and so with small 
corrections can be applied to the problem of a BEC confined in a ID harmonic trap 
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(Saito & Ueda 2003; Abdullaev et. al. 2003; Montesinos et. al. 2004). Note that for 
f{t) = one can derive a simple condition: J2VV - J3WW = J^V'"^ - J\JzW~'^ + 
4eF^ exp(— y^). This equation depends on time only implicitly, and is used when 
describing collapse scenarios (Section 6). 

In the corresponding Q2D treatment we drop the z dimension in equation (2.1). 
Wc assume that in this direction the profile of the wavcfunction is fixed and re- 
produces the ground state ^0 of the single lattice cell or harmonic potential, as 
in Saito & Ueda (2003) and Montesinos et. al. (2004). The reduced potential in 
2D GP will take the form U{g,t) = f {t) {1/2) loIq^ . In the VA wc take V con- 
stant, Vb as found from the equation 4£/V^''exp (— V^^) = 1, and only solve equa- 
tion (3.6). In numerical simulations we rescale the nonlinear coupling coefficient 
52D = 5 X (/ |V'o|^^od^;)/(/V'od2:). 

4. Numerical results 

We simulated both the full GPE, equation (2.1), using an axisymmetric code (for 
3D), a Cartesian code (for Q2D), and the variational equations for comparison. 
Numerical simulations followed the path outlined in figure 1. 

An example of the numerical results in the moderate confinement regime is 
shown in figure 2. The parameters used in the simulations would correspond, for 
^^Rb atoms, to an OL period of A = 1.5 jim., an initial radial-confinement frequency 
of u)± = 2-iT X 160 Hz, an FRM frequency oi Q = 2-it x 19 kHz, a lattice depth of 
£/ = 25 recoil energies, an effective nonlinear coefficient of Na = ±2 • 10~^ m, 
where N is the number of atoms per lattice cell, with a total number of atoms 
in the range of 10'' — 10^. The corresponding values of the normalized parameters 
are given in the figure captions. This figure shows the evolution of the central- 
peak's amplitude versus time in dimensionless units, defined by equation (2.1). 
After an initial transient, a stable breathing structure is established. The difference 
between the dynamics in the 3D treatment (lower curve) and the corresponding Q2D 
treatment (upper curve) is obvious. This is an example of a moderate confinement, 
when only the 3D treatment describes the dynamics of the system properly. The 
inset shows the 3D soliton structure for a fixed moment of time. 

The dynamics in the weak confinement regime are displayed in figure 3. The 
figure reveals the infiuence of neighbouring solitons on the amplitude of the central 
peak. The thick curve corresponds to the evolution of the full multipeak structure, 
similar to that shown in the inset of figure 2. The interaction between neighbouring 
solitons can be seen in the form of a beat. This phenomenon can be explained by 
the fact that the oscillation periods of the adjacent solitons are slightly different as 
there is a small difference in the number of atoms in neighbouring cells. To obtain 
the thin curve we repeated the above calculations up to the time t = 7500 and 
then removed all but the central peak. The structure so formed can be seen the 
inset of figure 3. The beat is no longer visible. The difference between these two 
cases proves that the dynamics of stable breathers as described here is a collective 
multi-peak phenomenon. 

In figure 4 we collected results of a systematic scan of parameter space based 
on GPE simulations and compared them with the predictions of the VA (a similar 
analysis can be performed if we replace ID OL with a ID harmonic trap - the 
conclusion does not depend on the form of confinement in the z direction). The 
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Figure 2. Comparison of the central peak amplitude evolution in the 3D treatment (lower 
curve) and the corresponding Q2D treatment (upper curve). This is an example of a 
moderate confinement, when only a fully 3D treatment describes the dynamics of the 
system properly. The normalized parameters are gof = —22, gu = 4gof, ef = 50, O = 36, 
LJ± — 0.3, ti = 30, t2 — 100, ts — 120, and t4 — 130. The inset shows the structure in 
the 3D treatment (in the Q2D treatment z dimension is dropped). The unit of time is 
/{-n^h). Stabilization in the z direction appears to be by attractive interaction. The 
soliton in this direction is discrete rather than gap-type. 

0.4 1 1 1 . 1 1 1 



0.3- 




t 



Figure 3. The thick and thin lines show the evolution of the amphtude of the central peak, 
respectively, of the three-peaked soliton, and in the case when, at t = 7500, the side peaks 
were suddenly removed (the latter configuration is shown in the inset). This is an example 
of a weak confinement, when individual solitons do interact. The normalized parameters 
are got = -18, git = 4gof, ti = 20.5, Q, = 22, ui± = 0.3, ti = 30, t2 = 100, ta = 120, and 
t4 = 130. 

agreement between VA and direct simulations is very good. The borders of the VA 
stability regions were found analytically e. g. upon considering effective potentials. 
A detailed derivation is presented in §5. As seen from figure 4 a), in a fully 3D 
treatment we have found two islands of stability. Note the similarity of the lower 
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Figure 4. Stability regions for solitons in the (Iffo/I, fi) plane, as predicted by the variational 
approximation (shaded area) , and found from direct simulations of the Gross-Pitaevskii 
equation (circles). Both a) corresponding to €f = 50 and flo = 26.76 and c) corresponding 
to Ef — 200 and flo — 55.06 were obtained from a 3D analysis and b) is the result of a 
Q2D treatment. The other parameters are as in figure 2. The lower regions in a) and c) 
correspond to Q2D solitons and those corresponding to the upper regions are fully 3D. 
As we see, as the confinement increases, the lower region obtained from the 3D analysis 
becomes more and more like that following from Q2D. The curves limiting stability regions 
were obtained analytically, see equations (5.5), (5.6), (5.11), and (5.17). The light-gray area 
in a) and c) is the stability region in the Q2D treatment appropriately rescaled. 



regions of a) and c) to that of figure 4 b), which portrays the results of the Q2D 
treatment. This region corresponds to Q2D solitons. On the other hand, the upper 
region contains fully 3D solitons. They appear when the frequency of modulation 
exceeds the lowest excitation frequency of the confining potential. As wc sec in 
figure 4, when the strength of the lattice £/ is increased, this region moves towards 
higher frequencies (and has risen beyond the scope of c)), the Q2D region expands 
upwards, and the whole picture becomes more and more like figure 4 b). 
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5. Stability analysis 

We base our calculations on the variational equations (3.6), (3.7) with f{t) = 0. In 
the context of these equations, we will derive approximate formulas for borders of 
stability regions in parameter space. Several types of instability will be considered. 
Some of them can only be described by the full 3D treatment. All stability conditions 
derived below lead to limiting curves in figure 4, with surprising accuracy. They 
arc in good agreement with the results of numerical simulations of both variational 
equations (3.6), (3.7) and the GPE (2.1). 

Q2D treatment. In the Q2D treatment we neglect the dynamics of V, assuming 
V = Vo and only consider the variational equation for W, (3.6) 

W = ^[^A + Bcosint)]. (5.1) 

Here we introduced the constants A = Ji{go/gc—^) and B = —Jigi/gc- The critical 
nonlinearity gc = — Vb J1/J2 is equal to the threshold for collapse in the absence of 
rapid modulations. 

In the adiabatic approximation we separate W into slowly and rapidly oscil- 
lating parts W = w + S, and assume that the latter is relatively small, S <^ w. 
Equation (5.1) takes the form 

w + (^-^][-A + Bcos{nt)]. (5.2) 

In the adiabatic limit of small ^ = |ti;|/(wO) we obtain approximately 

^ = — s?^cos(m), 

A 35^ 

^ = -^^^^M-- (^-^^ 

The second formula is obtained using the first and can be interpreted as an equation 
of motion of a particle in an effective potential 

^ef = -^ + 4. (5-4) 

where C = Ji {go/gc — 1) /2 and D = [Jigi/{2gc^)]^ . 

We denote the initial width wo and choose wo = 0. From (5.4) we see that the 
particle will be trapped by the potential provided that f/ef (wo) < 0, as t/gf (00) = 0. 
This gives a stability condition 



Wo > l^-^j , (5.5) 

which leads to the lower limiting curve in figure 4 b). Note that the above requires 
C > 0, or equivalently go/gc > 1- This was the only stability condition known so 
far (Saito & Ueda 2003; Montesinos et. al. 2004). 

We now consider if the adiabatic assumption used when deriving equations (5.3) 
is consistent with the result, i.e. whether ^, defined above, remains small. It's max- 
imum value can easily be estimated. If the oscillations in the potential (5.4) are 
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small, then is small and so is ^. If on the other hand these oscillations are large, 
one cycle can be split into two stages, (a) motion in the potential C/ef ~ —C/w^ for 
large w, and (b) motion in the potential Ucf ~ Djufi for small w. The ^ coefficient 
reaches its maximum value during the (b) stage. All possible trajectories in the (b) 
stage potential are approximately described by a function, which scales with D and 
the initial conditions in the following way 

^max ~ \/^n"^w;;t„(M;o), (5.6) 

where w^in is the value of w at the smaller turning point and can be calculated from 
Wo, using equation (5.4). It turns out that if the right-hand side of (5.6) exceeds 
a critical value, an instability occurs in the numerical simulations. Wc observe an 
increase in the particle energy (equation (5.3) is no longer valid) around the smaller 
turning points, and consequently particle escape from the potential well (5.4). We 
stress that, for ^max less than critical, the system is stable for a very long time. 
Thus wc obtain the upper limit in figure 4 b). 

3D treatment. Now wc extend our considerations to the case when the axial 
width of the peak, V is no longer constant and this influences stabilization. We 
keep the assumption that Y is close to the width of the localized Wannier function 
of our lattice potential. Thus we write = Vq + ?7, where r; -C Vg. 

The value of Vb is taken such that the first two terms on the right hand side of 
equation (3.7) cancel (there are in fact two such values and we take the lower one, 
known to be stable). From the third, nonlinear term we take only the oscillating 
part (the first two non-oscillating linear terms are assumed to dominate). In first 
order in ry/Vo we end up with a harmonic oscillator type formula with a small 
driving term 

y = r? w -0^77 + 7 cos(m), (5.7) 

where ^ 'ilV^^^E{\-1Y^)eyi^{-V^) = {2/Vif-2/Vi and7 = hgi/iViW). 
The solution is 

'7=^^j^cos(l]t). (5.8) 

Now we rewrite equation (3.6), neglecting terms of higher order in tj/Vq 
^_ J2(ffo+gicos(m)) / 77 \ 

and substitute rj from (5.8). In the adiabatic limit, we repeat the reasoning that 
led to equation (5.4), to obtain an equation of motion for a particle in an effective 
potential, including a new term 

C/ef = --^ + 4 + 4' (5-10) 

w' w'^ 

where F = J2Jzgi/ \8Vq{Vl'^ - Og)]. If the last term in (5.10) can be neglected 
as compared to the middle one (this is true if O is close to Oq), another stability 
condition can be found 

1/2 



«;o > ( ^ ) . (5.11) 
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This gives the lower limit of the upper region in figure 4 a). Strictly speaking, 
all terms in (5.10) should be taken into account when calculating these stability 
conditions. However, approximate formulas (5.5) and (5.11) prove to bo surprisingly 
accurate for the range of parameters studied here. A condition analogous to (5.6) 
will follow: 

Ua. ~ \/FO-itz;-3„(«;o), (5.12) 

However, another, lower border of the stability region can be found when consid- 
ering resonance with double frequency of the forcing term. In this case we have; to 
include higher, anharmonic terms in equation (5.7). In the following simple reason- 
ing we will just take the first, r]^ term, 

jj « -nl'ij + aT]'^ +^cos{nt), (5.13) 

where a = 8Vq~^ - 2Vq~^. We define f2 = 20o - £• The linear solution of 5.13 is 

= -^cos[(21)o-e)t], (5.14) 

and the next order correction r]^^\ obtained when the quadratic term is included, 
satifies 

i)(i)+f72,7(i) =2ar;(")r;(i). (5.15) 
For a solution in the form r]^^^ = 6cos[(f2o — e/2)t] we get 

-(Oo-e/2)2 + 0^ = -||. (5.16) 

Thus the threshold for appearance of a strong resonance when the driving frequency 
n comes close to 2rio is given by 

where we take 7 = Jsgi/ (VqW^^^^) and the minimal value of the width during the 
evolution as Wmin = \/D/C. In fact Wmin cannot be any smaller, as it would imply 
another instability (compare to equation 5.5). A quadratic in Q is obtained. Simple 
manipulations show that e tends to zero as C tends to zero when gof = Qoc ~ 20 
and f2 = 2fio on figure 4 a). 

A fuller analysis would include an Tf" term in (5.13) and the result would depend 
on b, the amplitude of This is why our e (5.17) is just a threshold value. For 
a full analysis, see Landau & Lifshitz (1960). Thus we obtain the upper limiting 
curve figure 4 a). The 77^ term indicates that the resonance in question appears 
above this limiting curve in figure 4 a). 

Interestingly, the instability discussed in the Q2D analysis, associated with a 
breakdown of the adiabatic approximation (5.6), is strongly suppressed in the 3D 
treament in the range < < 2fio (see equation 5.12). This is due to flattening 
of the eff'ective potential (5.10) by the additional term, and a decrease of the value 
of the ^ coefficient during the evolution. 
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6. Collapse and spreadout 

When considering possible collapse scenarios we look at the leading terms in equa- 
tion (3.6), when r = tcoi — t tends to zero. The sin(r^i) term either tends to a 
constant or else to ztfir, in which latter case it can be ignored in our considera- 
tions. We consider the leading small r terms of equation (3.6) to be 

W= ^'^^:^J^\ T = -J2{go + gism{^t,oi)) (6.1) 
and of equation (3.7) 

V=^--^^^^, r' = -Js{go+gism{nt,oi)). (6.2) 

Now r and V can take either common sign, but we will see that only F, F' > can 

lead to collapse. As t ^ cither all three terms in equation (6.1) cancel, or else 
two predominate and cancel in lowest order r. When all three cancel, W ~ r^/^, 
V — > Vcoi (1)- When two terms cancel, either the first and third are of lowest order 
and W ^ V ^ t^I^ (2), or else the second and third cancel because Vcoi = ^ / Ji 
and we find that ~ y - ~ r^/^ (3). We now look at all three possibilities in 
some detail. 

(1) All three terms in equation (6.1) are of equal strength. We find from both 
equations (6.1) and (6.2) 

W = aT'/^ + -^T^/^[ilnrf-3lnr\, (6.3) 
V = Feoi + /3(t In r-r), (6.4) 

where a = V^iT/Vcoi - JiY'\ P = r'/ia^V^^), T > VcoiJi- 

(2) When the first and the third term in equation (6.1) cancel, we find 

W = UT^/^ +aT^/^ + 0{t^/^), (6.5) 
V = pT^/^ +pT^/^ +0{t'^/^), (6.6) 

where a = (25/6)i/5(r3/F')^/i°, /? = [25F'V(6F)]V5, r,F' > 0, and 

a = (5/28) (3/?-^ - SJio;-^) , (6.7) 
p = (5/28) (6JiQ!-^ - 11/?-^) . (6.8) 

(3) The second and third term in equation (6.1) cancel when V = T / Ji+ fir^/^ , 
W = aT2/3. However, we find that ^ -9Ji/3/(2r) and ^ gp' j2/(2/3r2). As 
Ji > 0, and r,r' take the same sign, no real a, /3 pair satisfies these equations. 
This possibility is ruled out on the level of the coefficients. 

On the other hand, spreadout, when it occurs, is such that W and V are pro- 
portional to t at large times. Alternatively, just W ^ t and Vq const. 

7. Conclusions 

Only for strong confinement is the Q2D analysis an adequate approach. For moder- 
ate confinement the 3D region of stability is richer from that of Q2D (see figure 4). 
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A large, new region of stability appears in the 3D treatment as compared to the 
Q2D treatment. This is proof of the fact that our solitons are truly 3D. We are 
more used to the effect of adding a new dimension simply shrinking or abolishing 
the basin of stability of solitons or waves. For example, water waves are unstable 
with respect to perturbations along their direction of propagation only when the 
depth exceeds a critical value (Benjamin & Fcir 1967). When, however, two dimen- 
sional perturbations are allowed, there will always be an unstable angle regardless 
of the depth (Hayes 1973). Another example is that of ID solitons of the Nonlinear 
Schrodingcr equation (NLS) with constant coefficients. They arc stable in ID, but 
unstable in 2D or 3D. This is also true for some NLS waves (Anderson et. al. 1979; 
Infeld & Rowlands 1980, 2000). In the problem treated here this is not the case. 
For situations corcsponding to much of the quasi two dimensional stability chart 
adding a degree of freedom stabilizes the soliton solution. The key to this dichotomy 
is the presence of a periodic modulation, absent in the above mentioned classical 
examples. This can be illustrated by a simple case involving oscillators. Take as the 
one dimensional version a forced oscillator problem: 

X -\- WqX = ycos{ijjQt). (7.1) 

If y is fixed, the solution has a secular component 2/t/(2a)o) sin(a)oi), and so the 
amplitude will grow as t. If however we allow a second degree of freedom, such that 
y also oscillates (for instance y + e^y = 0) the solution stabilizes, unless e = ±2a;o. 
In general this can be the case when there arc periodic modulations. This fact, 
obvious in oscillator theory, is perhaps less well known in the soliton context. 

The main result of this paper is pointing out the possibility of creating fully 3D 
breather solitons in a BEC confined by a ID optical lattice potential, corresponding 
to the upper region in figure 4. The stable patterns may feature a multi-cell struc- 
ture. Depending on the strength of the confinement we identified three different 
cases: a) strong confinement; practically no evolution in the z direction, thus Q2D 
can be applied, b) moderate confinement; fully 3D soliton-train, but no interaction 
between cells, c) weak confinement; solution in form a set of weakly interacting fun- 
damental solitons. The scheme proposed here is based on a combination of FRM 
and a ID optical lattice, and could be implemented in an experiment relatively 
easily. This would open the way to the creation of robust 3D solitons (breathers) 
in BECs. 
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